Â
  Â
In Part 2 of the AusCover section of the workshop, we will investigate the Changes in Normalized Difference Vegetation Index (NDVI). In particular we will cover the following points:
For conciseness the messages returned by R are not displayed (i.e. only the results).
The images and plots in this report can appear small at times. However, they look fine on a computer monitor when the code is run in a computer.
In this section we will explore the change in NDVI in an area west of Tara where Coal Seam Gas development has been happening in recent years. NDVI can be used to estimate the density of green vegetation on a patch of land. It is derived from remote sensing data, typically a space platform. Chlorophyll, the pigment in plant leaves, strongly absorbs red and blue light. Leaves, on the other hand, strongly reflect near-infrared light (NIR) and green light (this is why vegetation looks green to our eyes). NDVI takes advantage of this difference in behaviour for different light wavelengths by green vegetation, combining measurements at different wavelengths in this formula: NDVI = (NIR-Red)/(NIR+Red). NDVI ranges from -1 to 1. With different values indicating different situations:
Generally, to investigate NDVI change over time an atmospheric correction needs to be performed. The data that we use in this workshop has already been corrected to reduce contamination by cloud and other problems.
We start by loading the required libraries, setting up a working directory, and cleaning up memory.
# Load Libraries
# ===============
library(dplyr)
library(reshape2)
library(stringr)
library(RColorBrewer)
library(ggplot2)
library(gridExtra)
library(RStoolbox)
library(sp)
library(rgdal)
library(raster)
library(rasterVis)
##library(RgoogleMaps)
library(ggmap)
#library(maps)
#library(mapdata)
#library(maptools)
# Set Data Path (TERN AusCover data repository)
# =============================================
data.path = "http://qld.auscover.org.au/public/data/landsat/surface_reflectance/aus"
# Optional Steps (remove comments to run them)
# ==============
# Setting up my Current Working Directory
# ---------------------------------------
#getwd()
#my.CWDir = "C:/Users/uqbblanc/Documents/TERN/CWDir"
#setwd(my.CWDir)
#getwd()
# Clean up Memory
# ---------------
#rm(list=ls())
#list.files()
We import and display a Stamen map to explore the area that we will be working with (West of Tara). The first attempt, using the extent of this area, leads to too detailed map. We add a plot where we zoom out by 2 degrees and now we can identify a nearby location (Dalby).
# Define the Extent of the Area to Explore
# ----------------------------------------
WofTara.extent = extent(150.3727, 150.4963, -27.0767, -26.9819)
WofTara.extent
## class : Extent
## xmin : 150.3727
## xmax : 150.4963
## ymin : -27.0767
## ymax : -26.9819
# Get Map with for West of Tara Extent
# ------------------------------------
WofTara.StamenMap = get_stamenmap(WofTara.extent[c(1,3,2,4)], maptype="terrain")
# Create Plot of the Map
# ----------------------
WofTara.Map.P1 =
ggmap(WofTara.StamenMap) + labs(title= "Study Area (W of Tara): Too detailed", x="Longitude", y="Latitude") +
theme(plot.title = element_text(hjust = 0.5, size=15, face="bold"),
axis.title = element_text(size=11, face="bold"), axis.text=element_text(size=9) )
#WofTara.Map.P1 # Plot individually for a bigger image
# Zoom out by 2 degrees & Create New Plot
# ---------------------------------------
# Download Map with new wider extent
WofTaraPlus2.StamenMap = get_stamenmap((WofTara.extent+2)[c(1,3,2,4)], maptype="terrain")
# Create Plot for new wider map
WofTaraPlus2.Map.P2 =
ggmap(WofTaraPlus2.StamenMap) + labs(title= "Study Area (W of Tara) + 2 Degrees", x="Longitude", y="Latitude") +
theme(plot.title = element_text(hjust = 0.5, size=15, face="bold"),
axis.title = element_text(size=11, face="bold"), axis.text=element_text(size=9) )
# Create Polygon with Original Extent and Add it to the Plot
WofTara.SP = as(WofTara.extent, "SpatialPolygons")
WofTara.SP # Has no CRS
## class : SpatialPolygons
## features : 1
## extent : 150.3727, 150.4963, -27.0767, -26.9819 (xmin, xmax, ymin, ymax)
## coord. ref. : NA
proj4string(WofTara.SP) = CRS("+init=epsg:4326") # Add CRS
WofTara.SP
## class : SpatialPolygons
## features : 1
## extent : 150.3727, 150.4963, -27.0767, -26.9819 (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
# Add Polygon to the Plot
WofTaraPlus2.Map.P2 = WofTaraPlus2.Map.P2 + geom_polygon(data=WofTara.SP, aes(x=long, y=lat), color="red", alpha=0)
#WofTaraPlus2.Map.P2 # Plot individually for a bigger image
# Show Both Plots
# ----------------
grid.arrange(WofTara.Map.P1, WofTaraPlus2.Map.P2, nrow=1)
This function includes many of the steps in the INTRODUCTORY STEPS section. Specifically the function perform these tasks:
The function takes the following arguments:
data.fpn: Data file path and namedl.fn: Name for the downloaded filerB.extent: RasterBrick extentrB.crs: RasterBrick coordinates reference systemThere are a number of print commands within the function ‘commented out’. These commands are not required for the correct functioning of the function. However, they have been included because they might make it easier to understand what the function is doing. To investigate further what the function is doing at each step simply remove the comment (#) symbol in front of the relevant print command.
get.B3B4B5.rB.f = function(data.fpn, dl.fn, rB.extent, rB.crs){
# Download the data file
# ======================
download.file(url=data.fpn, destfile=dl.fn, method='auto')
#print(list.files())
# Load the data file required bands
# =================================
LSSR.Band3.rL = raster(dl.fn, band=3)
#print(LSSR.Band3.rL)
LSSR.Band4.rL = raster(dl.fn, band=4)
#print(LSSR.Band4.rL)
LSSR.Band5.rL = raster(dl.fn, band=5)
#print(LSSR.Band5.rL)
# Subset Rasters
# ==============
# Project Extent to a CRS=EPSG:3577
# ---------------------------------
#print(rB.extent)
# Create a Polygon from Extent, Reproject the Poloygon and use its extension
rB.extent.GeoCoord.SP = as(rB.extent, "SpatialPolygons")
#print(rB.extent.GeoCoord.SP) # Has no CRS
proj4string(rB.extent.GeoCoord.SP) = CRS("+init=epsg:4326")
#print(rB.extent.GeoCoord.SP)
rB.extent.EPSG3577.SP = spTransform(rB.extent.GeoCoord.SP, CRS("+init=epsg:3577"))
#print(rB.extent.EPSG3577.SP)
rB.extent.EPSG3577 = extent(rB.extent.EPSG3577.SP)
#print(rB.extent.EPSG3577)
# Crop (i.e. subset) the raster layers to our area of interest
# ---------------------------------
LSSR.Subset.Band3.rL = crop(LSSR.Band3.rL, rB.extent.EPSG3577)
#print(LSSR.Subset.Band3.rL)
LSSR.Subset.Band4.rL = crop(LSSR.Band4.rL, rB.extent.EPSG3577)
#print(LSSR.Subset.Band4.rL)
LSSR.Subset.Band5.rL = crop(LSSR.Band5.rL, rB.extent.EPSG3577)
#print(LSSR.Subset.Band5.rL)
# Replace values < 0 with NAs
# ===========================
# Band 3: Red
# ------------
LSSR.Subset.Band3.rL[LSSR.Subset.Band3.rL < 0] = NA
#print(LSSR.Subset.Band3.rL)
#print(freq(is.na(LSSR.Subset.Band3.rL)))
#print(freq(is.na(LSSR.Subset.Band3.rL))/ncell(LSSR.Subset.Band3.rL))
# Band 4: NIR
# ------------
LSSR.Subset.Band4.rL[LSSR.Subset.Band4.rL < 0] = NA
#print(LSSR.Subset.Band4.rL)
#print(freq(is.na(LSSR.Subset.Band4.rL)))
#print(freq(is.na(LSSR.Subset.Band4.rL))/ncell(LSSR.Subset.Band4.rL))
# Band 5: SWIR
# ------------
LSSR.Subset.Band5.rL[LSSR.Subset.Band5.rL < 0] = NA
#print(LSSR.Subset.Band5.rL)
#print(freq(is.na(LSSR.Subset.Band5.rL)))
#print(freq(is.na(LSSR.Subset.Band5.rL))/ncell(LSSR.Subset.Band5.rL))
# Create a multi-layered RasterBrick object
# =========================================
LSSR.Subset.Bands3to5.rB = brick(LSSR.Subset.Band3.rL, LSSR.Subset.Band4.rL, LSSR.Subset.Band5.rL)
names(LSSR.Subset.Bands3to5.rB) = c("Band3.Red", "Band4.NIR", "Band5.SWIR")
#print(LSSR.Subset.Bands3to5.rB)
# Return RasterBrick with the Required Bands Subset to the desired Extent
# =======================================================================
return(LSSR.Subset.Bands3to5.rB)
} # get.B3B4B5.rB.f = function(data.fpn, dl.fn, rB.extent, rB.crs){
get.B3B4B5.rB.f) & Plot themUse our new function to get the Landsat Surface Reflectance datasets for Autumn 2007 and Autumn 2017, then subset the required bands and finally create and return a RasterBrick with the resulting RasterLayers. The RasterBricks for these periods are then plotted together.
getwd()
## [1] "C:/Users/uqbblanc/Documents/TERN/CWDir/FinalVersions_ToSend/SuperFinal/AusCover"
# Data Path and Data (Landsat Surface Refectance, LSSR) Files Names
data.path = "http://qld.auscover.org.au/public/data/landsat/surface_reflectance/aus"
data.2007.fn = "lztmre_aus_m200703200705_dbia2.vrt"
data.2017.fn = "l8olre_aus_m201703201705_dbia2.vrt"
# Create RasterBricks with LSSR Data for the required Extent and Period
# =====================================================================
# Autumn 2007
# -----------
# Create RasterBrick with data
wofTara.2007.LSSR.rB = get.B3B4B5.rB.f( data.fpn=paste(data.path, data.2007.fn, sep="/"), dl.fn="LSSR_2007.vrt", rB.extent=WofTara.extent, rB.crs=CRS("+init=epsg:4326") )
wofTara.2007.LSSR.rB
## class : RasterBrick
## dimensions : 408, 451, 184008, 3 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 1793835, 1807365, -3066365, -3054125 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
## data source : in memory
## names : Band3.Red, Band4.NIR, Band5.SWIR
## min values : 270, 938, 432
## max values : 3404, 4203, 5016
# Create Plot of RasterBrick
wofTara.2007.LSSR.p = ggRGB(wofTara.2007.LSSR.rB, r=3, g=2, b=1) +
labs(title= "LSSR: Autumn 2007", x="Longitude", y="Latitude") +
theme(plot.title = element_text(hjust = 0.5, size=20, face="bold"),
axis.title = element_text(size=12, face="bold"), axis.text=element_text(size=9) )
# Autumn 2017
# -----------
# Create RasterBrick with data
wofTara.2017.LSSR.rB = get.B3B4B5.rB.f( data.fpn=paste(data.path, data.2017.fn, sep="/"), dl.fn="LSSR_2017.vrt",
rB.extent=WofTara.extent, rB.crs=CRS("+init=epsg:4326") )
wofTara.2007.LSSR.rB
## class : RasterBrick
## dimensions : 408, 451, 184008, 3 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 1793835, 1807365, -3066365, -3054125 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
## data source : in memory
## names : Band3.Red, Band4.NIR, Band5.SWIR
## min values : 270, 938, 432
## max values : 3404, 4203, 5016
# Create Plot of RasterBrick
wofTara.2017.LSSR.p = ggRGB(wofTara.2017.LSSR.rB, r=3, g=2, b=1) +
labs(title= "LSSR: Autumn 2017", x="Longitude", y="Latitude") +
theme(plot.title = element_text(hjust = 0.5, size=20, face="bold"),
axis.title = element_text(size=12, face="bold"), axis.text=element_text(size=9) )
# Display both RasterBrick Plots in one plot
# ==========================================
grid.arrange(wofTara.2007.LSSR.p, wofTara.2017.LSSR.p, nrow=1)
In R raster calculations can be performed using:
'Higher Level Functions': First we place the formula with the required raster calculations in a function and then we use a Higher Level Function to invoke this function on relevant raster objects. This is more efficient (faster) option, which should be used for larger rasters and/or complex formulas. The Higher Level Function used depends on the number of raster object involved in the calculations:
calc and overlay are most commonly used.overlay is typically used.See TERN’s DSDP Tutorial ‘Using Raster Data in R’ for details.
To compute NDVI values we first create a formula with the required maths (for details on the NDVI see text at the beginning of this section (NORMALIZED DIFFERENCE VEGETATION INDEX (NDVI) CHANGE ANALYSIS). We actually create two formulas, the first one passing the RasterBrick to be used with calc and the second one passing individual layers to be used with overlay. Then we calculate this index in different ways:
calc and (3b) then using overlay.Sometimes calc and overlay don’t work; we show a trick (generating an empty raster with the right features and initialising it before assigning the values we need for our work to it) to help it work for us. The we compare the results of these 4 approaches to prove that the all produce the same results.
NOTE that in these datasets the value 32767 is used as a no-data code, so we will replace this value by NA
#=========================================================================================
# Create functions to Calculate NDVI
#=========================================================================================
# Function 1: Passing the RasterBrick
# -----------------------------------
calc.NDVI.f1 = function(rb) {
NDVI.rL = (rb[["Band4.NIR"]] - rb[["Band3.Red"]]) / (rb[["Band4.NIR"]] + rb[["Band3.Red"]])
return(NDVI.rL)
} # calc.NDVI.f1 = function(rb) {
# Function 2: Passing both RaterLayers
# -----------------------------------
calc.NDVI.f2 = function(rl.B4, rl.B3) {
NDVI.rL = (rl.B4 - rl.B3) / (rl.B4 + rl.B3)
return(NDVI.rL)
} # calc.NDVI.f2 = function(rl.B4, rl.B3) {
#=========================================================================================
# (1) Calculate NDVIs using Raster Algebra (= Raster Maths)
#=========================================================================================
# NDVI for 2007
# -------------
wofTara.2007.NDVI.rL.1 = (wofTara.2007.LSSR.rB[["Band4.NIR"]] - wofTara.2007.LSSR.rB[["Band3.Red"]]) /
(wofTara.2007.LSSR.rB[["Band4.NIR"]] + wofTara.2007.LSSR.rB[["Band3.Red"]])
wofTara.2007.NDVI.rL.1
## class : RasterLayer
## dimensions : 408, 451, 184008 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 1793835, 1807365, -3066365, -3054125 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
## data source : in memory
## names : layer
## values : -0.08312343, 0.7602862 (min, max)
# NDVI for 2017
# -------------
wofTara.2017.NDVI.rL.1 = (wofTara.2017.LSSR.rB[["Band4.NIR"]] - wofTara.2017.LSSR.rB[["Band3.Red"]]) /
(wofTara.2017.LSSR.rB[["Band4.NIR"]] + wofTara.2017.LSSR.rB[["Band3.Red"]])
wofTara.2017.NDVI.rL.1
## class : RasterLayer
## dimensions : 408, 451, 184008 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 1793835, 1807365, -3066365, -3054125 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
## data source : in memory
## names : layer
## values : -0.6510903, 0.8544872 (min, max)
(wofTara.2017.LSSR.rB[["Band3.Red"]] - wofTara.2017.LSSR.rB[["Band4.NIR"]]) /
(wofTara.2017.LSSR.rB[["Band4.NIR"]] + wofTara.2017.LSSR.rB[["Band3.Red"]])
## class : RasterLayer
## dimensions : 408, 451, 184008 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 1793835, 1807365, -3066365, -3054125 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
## data source : in memory
## names : layer
## values : -0.8544872, 0.6510903 (min, max)
#=========================================================================================
# (2) Directly calling one of the functions we created above to compute NDVI
#=========================================================================================
# NDVI for 2007
# -------------
wofTara.2007.NDVI.rL.2 = calc.NDVI.f1(rb = wofTara.2007.LSSR.rB)
wofTara.2007.NDVI.rL.2
## class : RasterLayer
## dimensions : 408, 451, 184008 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 1793835, 1807365, -3066365, -3054125 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
## data source : in memory
## names : layer
## values : -0.08312343, 0.7602862 (min, max)
# NDVI for 2017
# -------------
wofTara.2017.NDVI.rL.2 = calc.NDVI.f1(rb = wofTara.2017.LSSR.rB)
wofTara.2017.NDVI.rL.2
## class : RasterLayer
## dimensions : 408, 451, 184008 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 1793835, 1807365, -3066365, -3054125 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
## data source : in memory
## names : layer
## values : -0.6510903, 0.8544872 (min, max)
#=========================================================================================
# (3) Using Higher Level Functions
#=========================================================================================
#-----------------------------------------------------------------------------------------
# (3a) Using 'calc' & our 1st function to calculate NDVI
#-----------------------------------------------------------------------------------------
# Section before 'TRICK' commented out because it doesnt' work. Can give it a try
# by removing the comments.
## NDVI for 2007
## .............
#wofTara.2007.NDVI.rL.3a1 = calc(rb=wofTara.2007.LSSR.rB, fun=calc.NDVI.f1)
#wofTara.2007.NDVI.rL.3a1
# NDVI for 2017
# .............
#wofTara.2017.NDVI.rL.3a1 = calc(rb=wofTara.2017.LSSR.rB, fun=calc.NDVI.f1)
#wofTara.2017.NDVI.rL.3a1
# TRICK: Create & Initialise the Raster container and then copy the RasterLayers on it
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# All years have the same dimensions
rL = raster( ncol=ncol(wofTara.2007.LSSR.rB[["Band4.NIR"]]),
nrow=nrow(wofTara.2007.LSSR.rB[["Band4.NIR"]]) )
rL.B4.NIR = init(rL, fun=runif)
rL.B3.Red = init(rL, fun=runif)
# NDVI for 2007
# .............
rL.B4.NIR = wofTara.2007.LSSR.rB[["Band4.NIR"]]
rL.B3.Red = wofTara.2007.LSSR.rB[["Band3.Red"]]
rB = brick(rL.B4.NIR, rL.B3.Red)
names(rB) = c("Band4.NIR", "Band3.Red")
wofTara.2007.NDVI.rL.3a2 = calc(rB, fun=calc.NDVI.f1)
wofTara.2007.NDVI.rL.3a2
## class : RasterLayer
## dimensions : 408, 451, 184008 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 1793835, 1807365, -3066365, -3054125 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
## data source : in memory
## names : layer
## values : -0.08312343, 0.7602862 (min, max)
# NDVI for 2017
# .............
rL.B4.NIR = wofTara.2017.LSSR.rB[["Band4.NIR"]]
rL.B3.Red = wofTara.2017.LSSR.rB[["Band3.Red"]]
rB = brick(rL.B4.NIR, rL.B3.Red)
names(rB) = c("Band4.NIR", "Band3.Red")
wofTara.2017.NDVI.rL.3a2 = calc(rB, fun=calc.NDVI.f1)
wofTara.2017.NDVI.rL.3a2
## class : RasterLayer
## dimensions : 408, 451, 184008 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 1793835, 1807365, -3066365, -3054125 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
## data source : in memory
## names : layer
## values : -0.6510903, 0.8544872 (min, max)
#-----------------------------------------------------------------------------------------
# (3b) Using 'overlay' & our 2nd function to calculate NDVI
#-----------------------------------------------------------------------------------------
# Section before 'TRICK' commented out because it doesnt' work. Can give it a try
# by removing the comments.
## NDVI for 2007
## .............
#wofTara.2007.NDVI.rL.3b1 = overlay( rl.B4=wofTara.2007.LSSR.rB[["Band4.NIR"]],
# rl.B3=wofTara.2007.LSSR.rB[["Band3.Red"]],
# fun=calc.NDVI.f2)
#wofTara.2007.NDVI.rL.3b1
## NDVI for 2017
## .............
#wofTara.2017.NDVI.rL.3b1 = overlay( rl.B4=wofTara.2017.LSSR.rB[["Band4.NIR"]],
# rl.B3=wofTara.2017.LSSR.rB[["Band3.Red"]],
# fun=calc.NDVI.f2)
#wofTara.2017.NDVI.rL.3b1
# TRICK: Create & Initialise the Raster container and then copy the RasterLayers on it
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# All years have the same dimensions
rL = raster( ncol=ncol(wofTara.2007.LSSR.rB[["Band4.NIR"]]),
nrow=nrow(wofTara.2007.LSSR.rB[["Band4.NIR"]]) )
rL.B4.NIR = init(rL, fun=runif)
rL.B3.Red = init(rL, fun=runif)
# NDVI for 2007
# .............
rL.B4.NIR = wofTara.2007.LSSR.rB[["Band4.NIR"]]
rL.B3.Red = wofTara.2007.LSSR.rB[["Band3.Red"]]
wofTara.2007.NDVI.rL.3b2 = overlay(rL.B4.NIR, rL.B3.Red, fun=calc.NDVI.f2)
wofTara.2007.NDVI.rL.3b2
## class : RasterLayer
## dimensions : 408, 451, 184008 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 1793835, 1807365, -3066365, -3054125 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
## data source : in memory
## names : layer
## values : -0.08312343, 0.7602862 (min, max)
# NDVI for 2017
# .............
rL.B4.NIR = wofTara.2017.LSSR.rB[["Band4.NIR"]]
rL.B3.Red = wofTara.2017.LSSR.rB[["Band3.Red"]]
wofTara.2017.NDVI.rL.3b2 = overlay(rL.B4.NIR, rL.B3.Red, fun=calc.NDVI.f2)
wofTara.2017.NDVI.rL.3b2
## class : RasterLayer
## dimensions : 408, 451, 184008 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 1793835, 1807365, -3066365, -3054125 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
## data source : in memory
## names : layer
## values : -0.6510903, 0.8544872 (min, max)
#=========================================================================================
# Compare Results (those we can compare)
#=========================================================================================
# NDVI for 2007
# .............
all.equal(wofTara.2007.NDVI.rL.1, wofTara.2007.NDVI.rL.2)
## [1] TRUE
all.equal(wofTara.2007.NDVI.rL.1, wofTara.2007.NDVI.rL.3a2)
## [1] TRUE
all.equal(wofTara.2007.NDVI.rL.1, wofTara.2007.NDVI.rL.3b2)
## [1] TRUE
# NDVI for 2017
# .............
all.equal(wofTara.2017.NDVI.rL.1, wofTara.2017.NDVI.rL.2)
## [1] TRUE
all.equal(wofTara.2017.NDVI.rL.1, wofTara.2017.NDVI.rL.3a2)
## [1] TRUE
all.equal(wofTara.2017.NDVI.rL.1, wofTara.2017.NDVI.rL.3b2)
## [1] TRUE
# Remove unneeded Raster Objects
# ------------------------------
ls(pattern="rL.2|rL.3")
## [1] "wofTara.2007.NDVI.rL.2" "wofTara.2007.NDVI.rL.3a2"
## [3] "wofTara.2007.NDVI.rL.3b2" "wofTara.2017.NDVI.rL.2"
## [5] "wofTara.2017.NDVI.rL.3a2" "wofTara.2017.NDVI.rL.3b2"
rm(list=ls(pattern="rL.2|rL.3"))
We now plot the NDVI rasters using the levelplot function from the rasterVis package, which produces nice plots with more information than the standard plot function in the raster package`. We plot low NDVI values in red, intermediate values in yellow, and high values in green (the higher NDVI the darker the green).
To plot the NDVI’s of both year we first create a RasterBrick with both RasterLayers. We could plot both RasterLayers side to side but by default this would produce smaller plots as the plots would include density plots on their margins, which are not strictly comparable. They would also include two scales; plotting a RasterBrick a single scale is shown.
#=========================================================================================
# Plot NDVI rasters
#=========================================================================================
# Create the color palette
# ------------------------
# Red: low NDVIs, Yellow: middle NDVIs; Green: high NDVIs (Darker Greens indicate Higher NDVIs)
NDVI.breaks = seq(0,1,by=0.01)
NDVI.cols = colorRampPalette(c("red", "yellow", "darkgreen"))(length(NDVI.breaks)-1)
# Create side by side plots
# -------------------------
wofTara.NDVIs.rB = brick(wofTara.2007.NDVI.rL.1, wofTara.2017.NDVI.rL.1)
names(wofTara.NDVIs.rB) = c("NDVI.Autumn 2007", "NDVI.Autumn 2017")
levelplot(wofTara.NDVIs.rB, at=NDVI.breaks, col.regions=NDVI.cols, main="NDVI", names.attr=c("Autumn 2007","Autumn 2017"))
Now we explore the change in NDVI in the 10 year period (2007 to 2017). First, we compute the change in NDVI and then we visualise it.
To explore the change in NDVI we substract the NDVI values of 2007 from 2017. Therefore, positive values indicate an increase in NDVI, 0 no change, and negative values a decrease in NDVI.
Differences in NDVI can be due to many factors. Typically, the largest difference is due to climatic effects. Here the overall larger NDVIs in 2017 than in 2007 are due to a much wetter start of the year in 2017 than in in 2007. However, there are also specific areas of change. There are more watering points (indicated by very low NDVI values) and there is more coal seam gas infrastructure.
To make changes in NDVI slightly easier to interpret we normalise NDVI values (i.e. (NDVI - mean(NDVI)) / std(NDVI) ), so that 0 represent the mean NDVI and changes are in standard deviations of NDVI. To compute the summary statics values of a raster object (e.g. mean and standard deviation in our case) we need to use the cellStats function in the raster package. See TERN’s DSDP Tutorial ‘Using Raster Data in R’ for details.
As when we calculated NDVI, NDVI changes (both in absolute NDIV values and in stdev. units) can be computed using:
calc and overlay.For normalising the NDVI values in a raster we can use the function normImage in the package RStoolbox. We could incorporate this function both in ‘Raster Algebra’ formula/function and in the function called by the ‘Higher Level Function’. Here we will use ‘Raster Algebra’ both specifying the complete formula and taking advantage of the function normImage.
# Difference in 'raw' NDVI values
# ===============================
wofTara.NDVI.Diff.2017to2007.rL = wofTara.2017.NDVI.rL.1 - wofTara.2007.NDVI.rL.1
wofTara.NDVI.Diff.2017to2007.rL
## class : RasterLayer
## dimensions : 408, 451, 184008 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 1793835, 1807365, -3066365, -3054125 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
## data source : in memory
## names : layer
## values : -1.229459, 0.606049 (min, max)
# Difference in normalise NDVI values
# ===================================
# Method 1: Specifiying the whole formula
# ---------------------------------------
wofTara.NDVI.NormDiff.2017to2007.rL.m1 =
((wofTara.2017.NDVI.rL.1 - cellStats(wofTara.2017.NDVI.rL.1,mean)) / cellStats(wofTara.2017.NDVI.rL.1,sd)) -
((wofTara.2007.NDVI.rL.1 - cellStats(wofTara.2007.NDVI.rL.1,mean)) / cellStats(wofTara.2007.NDVI.rL.1,sd))
wofTara.NDVI.NormDiff.2017to2007.rL.m1
## class : RasterLayer
## dimensions : 408, 451, 184008 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 1793835, 1807365, -3066365, -3054125 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
## data source : in memory
## names : layer
## values : -10.9933, 3.723755 (min, max)
# Method 2: Using the 'normImage' function in the 'RStoolbox' package
# -------------------------------------------------------------------
wofTara.NDVI.NormDiff.2017to2007.rL.m2 = normImage(wofTara.2017.NDVI.rL.1) -
normImage(wofTara.2007.NDVI.rL.1)
wofTara.NDVI.NormDiff.2017to2007.rL.m2
## class : RasterLayer
## dimensions : 408, 451, 184008 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 1793835, 1807365, -3066365, -3054125 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
## data source : in memory
## names : layer
## values : -10.9933, 3.723755 (min, max)
# Compare results by both Normalisation Methods
# ---------------------------------------------
all.equal(wofTara.NDVI.NormDiff.2017to2007.rL.m1, wofTara.NDVI.NormDiff.2017to2007.rL.m2)
## [1] TRUE
To visualise the changes in ‘raw’ and normalised NDVI we will examine:
bwplot from the package rasterVis. This function creates violin plots, a boxplot with a kernel density plot.densityplot from the package rasterVis: Simple easy plots.ggplot from the package ggplot2: This looks prettier and is more flexible, but requires more work.We will start by creating a RasterBrick object to pass to many of the plotting functions (e.g. the functions in the package rasterVis).
# =========================
# Create RasterBrick object
# =========================
# 'Raw' NDVIs: Already done above (see 'Plot NDVI rasters' section)
# ...........
wofTara.RawNDVIs.rB = brick(wofTara.2007.NDVI.rL.1, wofTara.2017.NDVI.rL.1)
names(wofTara.RawNDVIs.rB) = c("Raw.NDVI.2007", "Raw.NDVI.2017")
# Normalised NDVIs
# ................
wofTara.2007.NormNDVI.rL.1 = normImage(wofTara.2007.NDVI.rL.1)
wofTara.2017.NormNDVI.rL.1 = normImage(wofTara.2017.NDVI.rL.1)
wofTara.NormNDVIs.rB = brick(wofTara.2007.NormNDVI.rL.1, wofTara.2017.NormNDVI.rL.1)
names(wofTara.NormNDVIs.rB) = c("Norm.NDVI.2007", "Norm.NDVI.2017")
# ========
# Boxplots
# ========
# Change in 'Raw' NDVI
# ....................
RawNDVI.IndvYrs.bwplot = bwplot(wofTara.RawNDVIs.rB, main="Raw NDVIs")
# Change in Normalised NDVI
# .........................
NormNDVI.IndvYrs.bwplot = bwplot(wofTara.NormNDVIs.rB, main="Normalised NDVIs")
# ========================
# Scatter plots with loess
# ========================
# Change in 'Raw' NDVI
# ....................
RawNDVI.IndvYrs.splot = splom(wofTara.RawNDVIs.rB, main="Raw NDVIs", plot.loess=TRUE, xlab='')
# Change in Normalised NDVI
# .........................
NormNDVI.IndvYrs.splot = splom(wofTara.NormNDVIs.rB, main="Normalised NDVIs", plot.loess=TRUE, xlab='')
# -------------------------------------
# Combine the 4 Plots in a single Graph
# -------------------------------------
grid.arrange( RawNDVI.IndvYrs.bwplot, NormNDVI.IndvYrs.bwplot,
RawNDVI.IndvYrs.splot, NormNDVI.IndvYrs.splot,
nrow=2, top="NDVI - Individuals Years" )
# =============
# Density Plots
# =============
# ----------------------------------------------------------
# Method 1: Using 'densityplot' from the package 'rasterVis`
# ----------------------------------------------------------
# Change in 'Raw' NDVI
# ....................
RawNDVI.IndvYrs.dp = densityplot(wofTara.RawNDVIs.rB, main="Individual Raw NDVIs")
RawNDVI.Diff.dp = densityplot(wofTara.NDVI.Diff.2017to2007.rL, main="Change in Raw NDVIs (2007-2017)")
# Change in Normalised NDVI
# .........................
NormNDVI.IndvYrs.dp = densityplot(wofTara.NormNDVIs.rB, main="Individual Normalised NDVIs")
NormNDVI.Diff.dp = densityplot(wofTara.NDVI.NormDiff.2017to2007.rL.m2, main="Change in Normalised NDVIs (2007-2017)")
# Combine 4 Plots in a Graph
# ..........................
grid.arrange( RawNDVI.IndvYrs.dp, RawNDVI.Diff.dp, NormNDVI.IndvYrs.dp, NormNDVI.Diff.dp,
nrow=2, top="NDVI: Individuals Years & Change" )
# ---------------------------------------------------
# Method 2: Using 'ggplot' from the package 'ggplot2`
# ---------------------------------------------------
# Create a DataFrame with the required values
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Create a DataFrame with the required values to pass to function 'ggplot' in the package 'ggplot2'
# Required values: NDVIs for 2007, 2017, and their difference in both scales: raw and normalised.
RawNDVI.2007 = values(wofTara.2007.NDVI.rL.1)
RawNDVI.2017 = values(wofTara.2017.NDVI.rL.1)
RawNDVI.Diff = values(wofTara.NDVI.Diff.2017to2007.rL)
NormNDVI.2007 = values(wofTara.2007.NormNDVI.rL.1)
NormNDVI.2017 = values(wofTara.2017.NormNDVI.rL.1)
NormNDVI.Diff = values(wofTara.NDVI.NormDiff.2017to2007.rL.m2)
NDVIs.df = data.frame( RawNDVI.2007, RawNDVI.2017, RawNDVI.Diff,
NormNDVI.2007, NormNDVI.2017, NormNDVI.Diff )
#summary(NDVIs.df)
# Change in 'Raw' NDVI
# ~~~~~~~~~~~~~~~~~~~~
# Density plot of Raw NDVI values for Individual Years:
# .....................................................
# Subset Dataset and convert wide format to long format
RawNDVI.IndvYrs.dflf = melt(NDVIs.df[c("RawNDVI.2007","RawNDVI.2017")])
RawNDVI.IndvYrs.dp2 = ggplot(RawNDVI.IndvYrs.dflf, aes(x=value, fill=variable)) +
geom_density(alpha=0.25)
# Density plot of the change in Raw NDVI values
# .............................................
RawNDVI.Diff.dp2 = ggplot(NDVIs.df, aes(x=RawNDVI.Diff)) +
geom_density(alpha=0.25,color="darkgreen", fill="lightgreen")
# Change in Normalised NDVI
# ~~~~~~~~~~~~~~~~~~~~~~~~~
# Density plot of Raw NDVI values for Individual Years:
# .....................................................
# Subset Dataset and convert wide format to long format
NormNDVI.IndvYrs.dflf = melt(NDVIs.df[c("NormNDVI.2007","NormNDVI.2017")])
NormNDVI.IndvYrs.dp2 = ggplot(NormNDVI.IndvYrs.dflf, aes(x=value, fill=variable)) +
geom_density(alpha=0.25)
# Density plot of the change in Raw NDVI values
# .............................................
NormNDVI.Diff.dp2 = ggplot(NDVIs.df, aes(x=NormNDVI.Diff)) +
geom_density(alpha=0.25,color="darkgreen", fill="lightgreen")
# Combine the 4 Plots in a single Graph
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
grid.arrange( RawNDVI.IndvYrs.dp2, RawNDVI.Diff.dp2, NormNDVI.IndvYrs.dp2, NormNDVI.Diff.dp2,
nrow=2, top="NDVI: Individuals Years & Change" )
# ============================================================
# Plot of the Raster Layer with the Normalised Change in NDVI
# ============================================================
# They have very differnt scales (raw NDVI and normalised NDVI values)
# -------------------
# Change in Raw NDVI
# -------------------
#wofTara.NDVI.Diff.2017to2007.rL
# Create the color palette
# ........................
abs.max = max(abs(values(wofTara.NDVI.Diff.2017to2007.rL)), na.rm=TRUE)
# Red: low NDVIs, Yellow: middle NDVIs; Green: high NDVIs (Darker Greens indicate Higher NDVIs)
RawNDVI.Diff.breaks = seq(-abs.max,abs.max, by=0.01)
RawNDVI.Diff.cols = colorRampPalette(c("red", "yellow", "darkgreen"))(length(RawNDVI.Diff.breaks)-1)
# Create Raster Plots with Density Plots
# ......................................
wofTara.RawNDVIDiff.p = levelplot( wofTara.NDVI.Diff.2017to2007.rL, at=RawNDVI.Diff.breaks,
col.regions=RawNDVI.Diff.cols, main="Raw NDVI Change (2017 - 2007)" )
# -------------------------
# Change in Normalised NDVI
# -------------------------
#wofTara.NDVI.NormDiff.2017to2007.rL.m2
# Create the color palette
# ........................
abs.max = max(abs(values(wofTara.NDVI.NormDiff.2017to2007.rL.m2)), na.rm=TRUE)
# Red: low NDVIs, Yellow: middle NDVIs; Green: high NDVIs (Darker Greens indicate Higher NDVIs)
NormNDVI.Diff.breaks = seq(-abs.max,abs.max, by=0.01)
NormNDVI.Diff.cols = colorRampPalette(c("red", "yellow", "darkgreen"))(length(NormNDVI.Diff.breaks)-1)
# Create Raster Plots with Density Plots
# ......................................
wofTara.NormNDVIDiff.p = levelplot( wofTara.NDVI.NormDiff.2017to2007.rL.m2, at=NormNDVI.Diff.breaks, col.regions=NormNDVI.Diff.cols, main="Normalised NDVI Change (2017 - 2007)")
# -----------------
# Plot Raster Plots
# -----------------
# `grid.arrange` provides multiple plots per page. They look good on screen,
# but might be too small in the Workshop/Tutorial report.
grid.arrange(wofTara.RawNDVIDiff.p, wofTara.NormNDVIDiff.p, nrow=1) # Too Small for Tut.
#wofTara.RawNDVIDiff.p # Used if 'grid.arrange' leads to plots too small
#wofTara.NormNDVIDiff.p # Used if 'grid.arrange' leads to plots too small
Aside from areas of where infrastructure has been put in and a small patch of clearing in the south-east corner, there is an area of noticeable greening up in the north-east.
To examine what is happening here overtime we can plot a time series of the green fraction. The green fraction is one of the bands of the fractional cover product. This is another AusCover product hosted by TERN.
Then we fit a trend line to the green data over the period of interest.
The server hosts time series stacks of cover fractions in virtual files that make pulling out time series a breeze. However, individual bands could also be iterated without too many problems.
To find our bearings we start by displaying the Greening Area. We do this in different ways: * Download and Display a map of the Greening Area from Stamen Maps. * Plot Satellite images of our study area (west of Tara) for years 2007 & 2017 with a (red) box indicating the Greening Area. * Plot the Normalised NDVI Change image with a (blue) box indicating the Greening Area.
The plots for the Satellite (Surface Reflectance) images are ggplot2 plots. The Normalised NDVI Change images is a lattice plot. The procedure to add a polygon box to ggplot2 & lattice plots differ. We plot both types of plots to show how to add a polygon in each case.
# Define the Extent of the Area to Explore
# ========================================
GreeningArea.extent = extent(150.4436, 150.4458, -27.0062, -27.0037)
GreeningArea.extent
## class : Extent
## xmin : 150.4436
## xmax : 150.4458
## ymin : -27.0062
## ymax : -27.0037
# Plot Web Map for the Greening Area
# ==================================
# Get Map
# -------
GreeningArea.StamenMap = get_stamenmap(GreeningArea.extent[c(1,3,2,4)], maptype="terrain")
# Create Plot of the Map
# ----------------------
# Stamen map of this area is not very detailed
GreeningArea.P1 =
ggmap(GreeningArea.StamenMap) + labs(title= "Greening Area", x="Longitude", y="Latitude") +
theme(plot.title = element_text(hjust = 0.5, size=15, face="bold"),
axis.title = element_text(size=11, face="bold"), axis.text=element_text(size=9) )
GreeningArea.P1
# Plot Greening Area over Satellite and NDVI Change Images
# ========================================================
# Create Polygon with the Greening Area Extent
# --------------------------------------------
GreeningArea.SP = as(GreeningArea.extent, "SpatialPolygons")
GreeningArea.SP # Has no CRS
## class : SpatialPolygons
## features : 1
## extent : 150.4436, 150.4458, -27.0062, -27.0037 (xmin, xmax, ymin, ymax)
## coord. ref. : NA
proj4string(GreeningArea.SP) = CRS("+init=epsg:4326") # Add CRS
GreeningArea.SP
## class : SpatialPolygons
## features : 1
## extent : 150.4436, 150.4458, -27.0062, -27.0037 (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
GreeningArea.EPSG3577.SP = spTransform(GreeningArea.SP, crs("+init=epsg:3577"))
GreeningArea.EPSG3577.SP
## class : SpatialPolygons
## features : 1
## extent : 1801860, 1802113, -3057837, -3057529 (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:3577 +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
# Plot Greening Area on Satellite Images
# --------------------------------------
# Plots 'wofTara.2007.LSSR.p' and 'wofTara.2017.LSSR.p' were created above.
# 2007 Image
# ~~~~~~~~~~
GreeningArea.2007.p = wofTara.2007.LSSR.p +
geom_polygon(data=GreeningArea.EPSG3577.SP, aes(x=long, y=lat),
color="red", alpha=0)
# 2007 Image
# ~~~~~~~~~~
GreeningArea.2017.p = wofTara.2017.LSSR.p +
geom_polygon(data=GreeningArea.EPSG3577.SP, aes(x=long, y=lat),
color="red", alpha=0)
# Plot Satellite Images from 2007 & 2017 with a Red Box for the Greening Area
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# `grid.arrange` provides multiple plots per page. They look good on screen,
# but might be too small in the Workshop/Tutorial report.
grid.arrange(GreeningArea.2007.p, GreeningArea.2017.p, nrow=1)
#GreeningArea.2007.p # Used if 'grid.arrange' leads to plots too small
#GreeningArea.2017.p # Used if 'grid.arrange' leads to plots too small
# Plot Greening Area on NDVI Change Image
# ---------------------------------------
# Parameters 'NormNDVI.Diff.breaks' & 'NormNDVI.Diff.cols' were created above
# (when genrating plot 'wofTara.NormNDVIDiff.p'.
levelplot( wofTara.NDVI.NormDiff.2017to2007.rL.m2, at=NormNDVI.Diff.breaks,
col.regions=NormNDVI.Diff.cols,
main="Normalised NDVI Change (2017 - 2007) with Greening Area in blue box") +
layer(sp.polygons(GreeningArea.EPSG3577.SP, col="blue", fill=NA))
We start by downloading the vrt file containing the green cover fraction. VRT is format (Virtual Format) driver for GDAL (Geospatial Data Abstraction Library) that allows building virtual GDAL datasets from other GDAL datasets. ‘.vrt’ files contain ‘VRT’ descriptions of datasets typically saved in XML format. For further details see ‘GDAL Virtual Format Tutorial’ by GDAL organization here
NOTE: That the name of the vrt file might have changed since this document was written, as files get updated (i.e. they grow) with new data.
# Data Path and Name
data.path = "http://qld.auscover.org.au/public/data/landsat/seasonal_fractional_cover/fractional_cover/aus"
#data.fn = "lztmre_aus_s198712201705_dima2_green.vrt"
#data.fn = "lztmre_aus_s198712201805_dima2_green.vrt"
data.fn = "lztmre_aus_s198712201808_dima2_green.vrt" # It's changed again!!!
# Download the data file (if it doesn't work try method='wget' or functions in Library 'RCurl')
download.file(url=paste(data.path, data.fn, sep="/"), destfile='GreenCoverFraction.vrt', method='auto')
list.files()
## [1] "ESA18_Workshop_AusCover_Overview.docx"
## [2] "ESA18_Workshop_AusCover_P1.html"
## [3] "ESA18_Workshop_AusCover_P1.R"
## [4] "ESA18_Workshop_AusCover_P1.Rmd"
## [5] "ESA18_Workshop_AusCover_P1_NoStamen.html"
## [6] "ESA18_Workshop_AusCover_P1_NoStamen.R"
## [7] "ESA18_Workshop_AusCover_P1_NoStamen.Rmd"
## [8] "ESA18_Workshop_AusCover_P2.html"
## [9] "ESA18_Workshop_AusCover_P2.R"
## [10] "ESA18_Workshop_AusCover_P2.Rmd"
## [11] "ESA18_Workshop_AusCover_P2_files"
## [12] "ESA18_Workshop_AusCover_P2_NoStamen.html"
## [13] "ESA18_Workshop_AusCover_P2_NoStamen.R"
## [14] "ESA18_Workshop_AusCover_P2_NoStamen.Rmd"
## [15] "ESA18_Workshop_AusCover_P2_v6f.html"
## [16] "ESA18_Workshop_AusCover_P2_v6f.R"
## [17] "ESA18_Workshop_AusCover_P2_v6f.Rmd"
## [18] "ESA18_Workshop_AusCover_ReadMe.docx"
## [19] "GreenCoverFraction.vrt"
## [20] "LSSR_2007.vrt"
## [21] "LSSR_2017.vrt"
The green cover fraction vrt file currently contains 122 bands/layers. Four new bands are created and added each year. Bands 77-120 correspond to the period of interest, 2007-2017 (one band has already been added for 2018). You can check this (and the file structure) opening it in text editor such as Notepad.
We create a new function to Load, Subset (i.e. Crop), Recode (the no-data values) and Combine the bands of interest. This function is partially similar to the function ‘get.B3B4B5.rB.f’ that we created above. Then we call this function to create a RasterBrick containing the Green Cover Fraction bands of interest (bands 77-120 for the period 2007-2017) subsetted to the extent of the Greening Area we want to explore.
#-----------------------------------------------------------------------------------------
# Create a function to Extract, Subset, Recode, and Combine the Bands into a Raster Brick
#-----------------------------------------------------------------------------------------
getprep.Bands.rB.f = function(dl.fn, bands.v, rB.extent, rB.crs){
# Create cookie cutter: Extent projected to CRS=EPSG:3577
# =======================================================
#print(rB.extent)
# Create a Polygon from Extent, Reproject the Poloygon and use its extension
rB.extent.GeoCoord.SP = as(rB.extent, "SpatialPolygons")
#print(rB.extent.GeoCoord.SP) # Has no CRS
proj4string(rB.extent.GeoCoord.SP) = CRS("+init=epsg:4326")
#print(rB.extent.GeoCoord.SP)
rB.extent.reprj.SP = spTransform(rB.extent.GeoCoord.SP, rB.crs)
#print(rB.extent.reprj.SP)
rB.extent.reprj = extent(rB.extent.reprj.SP)
#print(rB.extent.reprj)
# For each Band
# =============
for ( band.cnt in bands.v ) {
#print(band.cnt)
# Load the data to RasterLayer object
# -----------------------------------
GFC.rL = raster(dl.fn, band=band.cnt)
#print(GFC.rL)
# Subset (i.e. crop) raster layer to area of interest
# ---------------------------------------------------
GFC.Subset.rL = crop(GFC.rL, rB.extent.reprj)
is.na(GFC.Subset.rL)
#print(GFC.Subset.rL)
# Replace values < 0 with NAs
# ---------------------------
# GFC.Subset.rL[GFC.Subset.rL < 0] = NA
# Add Raster Layer to Raster Brick
# --------------------------------
if ( band.cnt == bands.v[1] )
GFC.Subset.rB = brick(GFC.Subset.rL)
else
GFC.Subset.rB = addLayer(GFC.Subset.rB, GFC.Subset.rL)
} # for ( band.cnt in bands.v ) {
# Re-conver to RasterBrick object (if prefered)
# =============================================
# 'addLayer' (used in the for loop) returns a RasterStack object, so we need to re-convert
# object to a RasterBrick (if this is the prefered object class).
GFC.Subset.rB = brick(GFC.Subset.rB)
# Add Layer Names to Raster Brick
# ===============================
# Create a vector of Raster Layer names (including the Band, Year, and Quarter)
band.part = paste("Band", bands.v, sep="")
initial.year = 1988
date.part = as.character(formatC(initial.year+bands.v/4-0.25,digits=2,format="f")) # '-0.25 'cos statrs in band1 not band0
band.names = paste(band.part, date.part, sep="_")
# Add the band names to the RasterBrick
names(GFC.Subset.rB) = band.names
#print(GFC.Subset.rB)
# Return RasterBrick with the Required Bands Subset to the desired Extent
# =======================================================================
return(GFC.Subset.rB)
} # getprep.Bands.rB.f = function(dl.fn, bands.v, rB.extent, rB.crs){
#----------------------------------------------------------------------------------------
# Call the function to Create a RasterBrick with processed Bands for period 2007 to 2017
#----------------------------------------------------------------------------------------
GreeningArea.rB =getprep.Bands.rB.f( dl.fn='GreenCoverFraction.vrt', bands.v=77:120,
rB.extent=GreeningArea.extent, rB.crs=CRS("+init=epsg:3577") )
GreeningArea.rB
## class : RasterBrick
## dimensions : 11, 9, 99, 44 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 1801845, 1802115, -3057845, -3057515 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
## data source : in memory
## names : Band77_2007.00, Band78_2007.25, Band79_2007.50, Band80_2007.75, Band81_2008.00, Band82_2008.25, Band83_2008.50, Band84_2008.75, Band85_2009.00, Band86_2009.25, Band87_2009.50, Band88_2009.75, Band89_2010.00, Band90_2010.25, Band91_2010.50, ...
## min values : 100, 100, 100, 116, 121, 112, 109, 112, 113, 116, 120, 110, 118, 126, 111, ...
## max values : 124, 128, 138, 149, 150, 142, 136, 128, 127, 135, 145, 127, 140, 153, 146, ...
First we compute a vector with the mean Green Cover Fraction (GCF) in the Greening Area (GA) for each raster layers. The raster layer values are in scale 0-255, so we re-scale the mean value to scale of 0-100 to represent the Green Cover Fraction in Percentage.
Then we create a vector with the Times (Years and Quarters) for each of these raster layers. The function that we used to create the Raster Brick with all the relevant data computed the date for each layer and assigned it as part of the corresponding band (i.e. layer) name. Here, we take advantage of this and extract the Times from the bands names. We could also have calculated it in a similar fashion as it was done within the function.
# Compute the Mean GCF in the Greening Area for each of the Bands
# ================================================================
GCF.GA.Means = cellStats(GreeningArea.rB*100/255, mean)
GCF.GA.Means
## Band77_2007.00 Band78_2007.25 Band79_2007.50 Band80_2007.75
## 42.07170 43.42246 45.26837 52.84611
## Band81_2008.00 Band82_2008.25 Band83_2008.50 Band84_2008.75
## 53.65815 48.69479 47.51040 46.94791
## Band85_2009.00 Band86_2009.25 Band87_2009.50 Band88_2009.75
## 46.27055 48.67895 52.75500 45.90216
## Band89_2010.00 Band90_2010.25 Band91_2010.50 Band92_2010.75
## 49.40780 55.65855 50.83779 NaN
## Band93_2011.00 Band94_2011.25 Band95_2011.50 Band96_2011.75
## 59.07704 60.67340 52.73123 52.49752
## Band97_2012.00 Band98_2012.25 Band99_2012.50 Band100_2012.75
## NaN 55.62370 54.23648 49.41573
## Band101_2013.00 Band102_2013.25 Band103_2013.50 Band104_2013.75
## NaN 57.22321 52.92533 50.03763
## Band105_2014.00 Band106_2014.25 Band107_2014.50 Band108_2014.75
## 51.25767 56.03882 54.03446 49.59794
## Band109_2015.00 Band110_2015.25 Band111_2015.50 Band112_2015.75
## 53.66211 54.51773 54.23252 51.46762
## Band113_2016.00 Band114_2016.25 Band115_2016.50 Band116_2016.75
## 54.96534 57.83720 56.28045 54.83462
## Band117_2017.00 Band118_2017.25 Band119_2017.50 Band120_2017.75
## 55.67835 58.90275 56.31214 52.38661
# Create a vector with the Times for each of Bands
# ================================================
#names(GCF.GA.Means)
GCF.GA.Times = as.numeric( sub(".*_","", names(GCF.GA.Means)) )
# OR: as.numeric( sapply(strsplit(names(GCF.GA.Means),split="[_]"),"[",2) )
GCF.GA.Times
## [1] 2007.00 2007.25 2007.50 2007.75 2008.00 2008.25 2008.50 2008.75
## [9] 2009.00 2009.25 2009.50 2009.75 2010.00 2010.25 2010.50 2010.75
## [17] 2011.00 2011.25 2011.50 2011.75 2012.00 2012.25 2012.50 2012.75
## [25] 2013.00 2013.25 2013.50 2013.75 2014.00 2014.25 2014.50 2014.75
## [33] 2015.00 2015.25 2015.50 2015.75 2016.00 2016.25 2016.50 2016.75
## [41] 2017.00 2017.25 2017.50 2017.75
To explore the Change in Green Cover Fraction in the Greening Area in the study period (2007-2017) we create two types of plots:
poly), and (3) a locally weighted regression (loess). For prettier and more flexible graphs we use the ggplot function from the ggplot2` package.bwplot from the package rasterVis to create this type of plot.# Scatter-plot (with lines) using different Smooths
# =================================================
# 'ggplot' requires a dataframe, so we first create a dataframe containing the
# Green Cover Fractions (%) and Times (Years & Quarters).
# Create dataframe
# ----------------
GCF.GA.df = data.frame(Time=GCF.GA.Times, Mean=GCF.GA.Means)
summary(GCF.GA.df)
## Time Mean
## Min. :2007 Min. :42.07
## 1st Qu.:2010 1st Qu.:49.42
## Median :2012 Median :52.85
## Mean :2012 Mean :52.35
## 3rd Qu.:2015 3rd Qu.:55.62
## Max. :2018 Max. :60.67
## NA's :3
# Create base plot (base + points)
# --------------------------------
base.p = ggplot(GCF.GA.df, aes(x=Time, y=Mean)) +
theme(plot.title=element_text(hjust = 0.5, size=14, face="bold")) +
labs(x="Time (Year.Quarter)", y="Green Cover Fraction (%)") +
scale_x_continuous(breaks=2007:2018) +
geom_line(color="blue")
#base.p
# Add different types of Smooths to the Base Plot
# -----------------------------------------------
# No Trend
base.NoSmooth.p = base.p + ggtitle("Green Cover Fraction in Greening Area (2007 to 2017)")
#base.NoSmooth.p
# Linear Model (LM) Fit
base.SmoothLM.p = base.p + ggtitle("Green Cover Fraction + LM Trend") +
stat_smooth(method = "lm", formula = y ~ x, size = 1)
#base.SmoothLM.p
# LM with 2nd Degree Polynomial Fit (using the function 'poly')
base.SmoothPoly2D.p = base.p +
ggtitle("Green Cover Fraction + LM with 2 Deg. Poly. Trend") +
stat_smooth(method = "lm", formula = y ~ poly(x,2), size = 1)
#base.SmoothPoly2D.p
# Logically Weighted Regression (Loess)
base.SmoothLoess.p = base.p + ggtitle("Green Cover Fraction + LOESS Trend") +
stat_smooth(method = "loess", formula = y ~ x, size = 1)
#base.SmoothLoess.p
# Display all Scatter-plots (with Lines) with different types of Smooths
# ----------------------------------------------------------------------
# `grid.arrange` provides multiple plots per page. They look good on screen,
# but might be too small in the Workshop/Tutorial report.
grid.arrange(base.NoSmooth.p, base.SmoothLM.p, base.SmoothPoly2D.p, base.SmoothLoess.p,
nrow=2)
#base.NoSmooth.p # Used if 'grid.arrange' leads to plots too small
#base.SmoothLM.p # Used if 'grid.arrange' leads to plots too small
#base.SmoothPoly2D.p # Used if 'grid.arrange' leads to plots too small
#base.SmoothLoess.p # Used if 'grid.arrange' leads to plots too small
# Violin Plots (Boxplot + Density Plots)
# ======================================
# Rather than using the default tick labels (too long including the Band and Time), we
# use tick labels including only the Time (Year & Quarter).
bwplot( GreeningArea.rB*100/255,
main="Green Cover Fraction in the Greening Area (2007 to 2017)",
xlab="Time (Year.Quarter)", ylab="Green Cover Fraction (%)",
scales=list(x=list(labels=format(GCF.GA.Times,nsmall=2),rot=45)) )